% Mass statistical test that compares 3 methods of template matching:
% - Kunchenko's polynomials
% - Signal cross-correlation
% - Sum of squared differences squares
%
% template and signal initializing
step = 0.01;
% template - canonacal gaussian
[xpsi psi] = template3(step);
% trust level to Kunchenko and Cross-correlation method
KunchenkoCorrelationTrustLevel = .4:.1:.9;
% SSD theshold level
SSDThresholdLevel = .05;
% maximum derivation in
trustedDerivation = 0.1;
% number of gaussians in experiments
gaussiansNumber = 10;
% S - number of statistic iterations (size of samlpe)
S = 100;

% ================= clean signal matching =================================

% effectivity results initialization

KunchenkoFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
CorrelationFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
SSDFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));

KunchenkoWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
CorrelationWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
SSDWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));

for i = 1: length(KunchenkoCorrelationTrustLevel)
    for count = 1:S
        % statistic experiment
        % signal generating - sum of gaussians that are generated randomly
        [x sign means fmeans] = signal8(step, gaussiansNumber);
        % getting maximums of clean signal
        [xSignMax SignMax] = maxEdges(x, sign, .8 * min(fmeans));
        
        
        % ===== Kunchenko's polynomials ======================================
        [~, ~, e]=KunchenkoMod(sign, psi, 6);
        scaledKunchenkoEffectogram = kunchenkoEffectogramScaling(x, e, xpsi);
        % getting maximums of thresholded effectogram that are over trust level
        [trustedMaximums ~] = maxEdges(x, scaledKunchenkoEffectogram, KunchenkoCorrelationTrustLevel(i));
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMaximums);
        KunchenkoFindedGaussiansNumber(i) = KunchenkoFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        KunchenkoWrongGaussiansNumber(i) = KunchenkoWrongGaussiansNumber(i)...
            + (length(trustedMaximums) - numberFindedGaussians);
        
        % ===== Signal cross-correlation ==================================
        r = correlation(sign, psi);
        scaledCorrelationEffectogram = correlationEffectogramScaling(x, r, xpsi);
        % getting maximums of thresholded effectogram
        [trustedMaximums ~] = maxEdges(x, scaledCorrelationEffectogram, KunchenkoCorrelationTrustLevel(i));
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMaximums);
        CorrelationFindedGaussiansNumber(i) = CorrelationFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        CorrelationWrongGaussiansNumber(i) = CorrelationWrongGaussiansNumber(i)...
            + (length(trustedMaximums) - numberFindedGaussians);
        
        % ===== Sum of squared differences squares =======================
        ssd = SSD(sign, psi);
        scaledSSDEffectogram = SSDEffectogramScaling(x, ssd, xpsi);
        % getting local minimums of SSD effectogram
        [trustedMinimums ~] = minEdges(x, scaledSSDEffectogram, SSDThresholdLevel);
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMinimums);
        SSDFindedGaussiansNumber(i) = SSDFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        SSDWrongGaussiansNumber(i) = SSDWrongGaussiansNumber(i)...
            + (length(trustedMinimums) - numberFindedGaussians);
    end;
end;

% getting average values
KunchenkoFindedGaussiansNumber = KunchenkoFindedGaussiansNumber / S;
CorrelationFindedGaussiansNumber = CorrelationFindedGaussiansNumber / S;
SSDFindedGaussiansNumber = SSDFindedGaussiansNumber / S;

KunchenkoWrongGaussiansNumber = KunchenkoWrongGaussiansNumber / S;
CorrelationWrongGaussiansNumber = CorrelationWrongGaussiansNumber / S;
SSDWrongGaussiansNumber = SSDWrongGaussiansNumber / S;

% plotting results
figure;
plot(KunchenkoCorrelationTrustLevel, [KunchenkoFindedGaussiansNumber; ...
    CorrelationFindedGaussiansNumber; ...
    SSDFindedGaussiansNumber]);
grid on;
title('Template Matching methods comparing for clean signal');
legend('Kunchenko''s polynomials', 'Signal cross-correlation', 'Sum of squared differences');
xlabel('trust level to Kunchenko and Cross-correlation method');
ylabel('Number of found gaussians');

figure;
plot(KunchenkoCorrelationTrustLevel, [KunchenkoWrongGaussiansNumber; ...
    CorrelationWrongGaussiansNumber; ...
    SSDWrongGaussiansNumber]);
grid on;
title('Template Matching methods comparing for clean signal');
legend('Kunchenko''s polynomials', 'Signal cross-correlation', 'Sum of squared differences');
xlabel('trust level to Kunchenko and Cross-correlation method');
ylabel('Number of wrong-found gaussians');

% ================= noised signal matching ================================
% noise level
SNR = 30;

% effectivity results initialization

KunchenkoFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
CorrelationFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
SSDFindedGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));

KunchenkoWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
CorrelationWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));
SSDWrongGaussiansNumber = zeros(1, length(KunchenkoCorrelationTrustLevel));

for i = 1: length(KunchenkoCorrelationTrustLevel)
    for count = 1:S
        % statistic experiment
        % signal generating - sum of gaussians that are generated randomly
        [x sign means fmeans] = signal8(step, gaussiansNumber);
        % getting maximums of clean signal
        [xSignMax SignMax] = maxEdges(x, sign, .8 * min(fmeans));
        % adding white noise to signal
        noised_sign = awgn(sign, SNR);
        
        % ===== Kunchenko's polynomials ======================================
        [~, ~, e]=KunchenkoMod(noised_sign, psi, 6);
        scaledKunchenkoEffectogram = kunchenkoEffectogramScaling(x, e, xpsi);
        % getting maximums of thresholded effectogram that are over trust level
        [trustedMaximums ~] = maxEdges(x, scaledKunchenkoEffectogram, KunchenkoCorrelationTrustLevel(i));
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMaximums);
        KunchenkoFindedGaussiansNumber(i) = KunchenkoFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        KunchenkoWrongGaussiansNumber(i) = KunchenkoWrongGaussiansNumber(i)...
            + (length(trustedMaximums) - numberFindedGaussians);
        
        % ===== Signal cross-correlation ==================================
        r = correlation(sign, psi);
        scaledCorrelationEffectogram = correlationEffectogramScaling(x, r, xpsi);
        % getting maximums of thresholded effectogram
        [trustedMaximums ~] = maxEdges(x, scaledCorrelationEffectogram, KunchenkoCorrelationTrustLevel(i));
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMaximums);
        CorrelationFindedGaussiansNumber(i) = CorrelationFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        CorrelationWrongGaussiansNumber(i) = CorrelationWrongGaussiansNumber(i)...
            + (length(trustedMaximums) - numberFindedGaussians);
        
        % ===== Sum of squared differences =======================
        ssd = SSD(sign, psi);
        scaledSSDEffectogram = SSDEffectogramScaling(x, ssd, xpsi);
        % getting local minimums of SSD effectogram
        [trustedMinimums ~] = minEdges(x, scaledSSDEffectogram, SSDThresholdLevel);
        % getting number of finded gaussians
        numberFindedGaussians = searchEffectivity(xSignMax, trustedDerivation, trustedMinimums);
        SSDFindedGaussiansNumber(i) = SSDFindedGaussiansNumber(i) + numberFindedGaussians;
        % getting number of wrong-found gaussians
        SSDWrongGaussiansNumber(i) = SSDWrongGaussiansNumber(i)...
            + (length(trustedMinimums) - numberFindedGaussians);
    end;
end;

% getting average values
KunchenkoFindedGaussiansNumber = KunchenkoFindedGaussiansNumber / S;
CorrelationFindedGaussiansNumber = CorrelationFindedGaussiansNumber / S;
SSDFindedGaussiansNumber = SSDFindedGaussiansNumber / S;

KunchenkoWrongGaussiansNumber = KunchenkoWrongGaussiansNumber / S;
CorrelationWrongGaussiansNumber = CorrelationWrongGaussiansNumber / S;
SSDWrongGaussiansNumber = SSDWrongGaussiansNumber / S;

% plotting results
figure;
plot(KunchenkoCorrelationTrustLevel, [KunchenkoFindedGaussiansNumber; ...
    CorrelationFindedGaussiansNumber; ...
    SSDFindedGaussiansNumber]);
grid on;
title(strcat('Template Matching methods comparing for noised signal SNR=', num2str(SNR), 'dB'));
legend('Kunchenko''s polynomials', 'Signal cross-correlation', 'Sum of squared differences');
xlabel('trust level to Kunchenko and Cross-correlation method');
ylabel('Number of founded gaussians');

figure;
plot(KunchenkoCorrelationTrustLevel, [KunchenkoWrongGaussiansNumber; ...
    CorrelationWrongGaussiansNumber; ...
    SSDWrongGaussiansNumber]);
grid on;
title(strcat('Template Matching methods comparing for noised signal SNR=', num2str(SNR), 'dB'));
legend('Kunchenko''s polynomials', 'Signal cross-correlation', 'Sum of squared differences');
xlabel('trust level to Kunchenko and Cross-correlation method');
ylabel('Number of wrong-found gaussians');